---
title: "Equity in NYC Park Access: A Spatial Analysis of Walking Distance"
subtitle: "STA 9750 Final Individual Report"
author: "Mirae Han"
date: "December 14, 2025"
format:
html:
toc: true
toc-depth: 3
toc-location: left
code-fold: true
code-summary: "Show code"
code-tools: true
number-sections: true
number-depth: 2
theme: cosmo
self-contained: true
fig-width: 10
fig-height: 6
df-print: paged
execute:
warning: false
message: false
cache: true
---
# Introduction {#sec-introduction}
## Group Overarching Question
**How equitable is access to public park space across different neighborhoods in New York City?**
Our team investigates park equity through three complementary dimensions:
1. **Quantity**: Park acreage per capita across neighborhoods
2. **Quality**: Park amenities and maintenance levels
3. **Proximity** (my focus): Walking distance to nearest park
## My Specific Question
**How does the average walking distance to the nearest park vary across the city and neighborhood demographic composition?**
### Research Focus
This analysis examines park accessibility through the lens of spatial proximity, measuring walking distances from residential areas to their nearest parks. By integrating demographic data (income and race/ethnicity), we assess whether park access is equitably distributed across NYC's diverse neighborhoods.
### Hypothesis
Based on environmental justice literature, we hypothesize that:
1. Lower-income neighborhoods have greater average distances to parks
2. Predominantly minority neighborhoods face reduced park accessibility
3. Spatial disparities exist across boroughs
---
# Data Acquisition and Processing {#sec-data}
```{r setup}
#| label: setup
#| include: true
# Clear environment
rm(list = ls())
invisible(gc())
# Required packages
required_packages <- c(
"tidyverse", # Data manipulation and visualization
"sf", # Spatial data handling
"viridis", # Color palettes
"scales", # Formatting
"jsonlite", # JSON parsing
"httr", # HTTP requests
"leaflet", # Interactive maps
"plotly", # Interactive plots
"patchwork" # Plot composition
)
# Function to install and load packages
install_and_load <- function(package) {
if (!require(package, character.only = TRUE, quietly = TRUE)) {
install.packages(package, dependencies = TRUE, repos = "https://cloud.r-project.org")
library(package, character.only = TRUE)
} else {
library(package, character.only = TRUE)
}
}
# Load all packages
invisible(sapply(required_packages, install_and_load))
# Setup data directory
data_dir <- "data/Group_Project"
if (!dir.exists(data_dir)) {
dir.create(data_dir, recursive = TRUE)
}
```
## NYC Parks Properties Data
**Source**: NYC Open Data - Parks Properties
**URL**: https://data.cityofnewyork.us/api/views/enfh-gkve/rows.csv
The NYC Parks Properties dataset contains over 1,800 park records with geographic boundaries as Well-Known Text (WKT) multipolygon strings. Key preprocessing steps included converting WKT to spatial features using `st_as_sf()`, mapping borough codes (M, B, Q, X, R) to full names, and filtering invalid geometries. All spatial data was transformed to EPSG:2263 (NYC State Plane, feet) for accurate distance calculations.
```{r load-parks}
#| label: load-parks
parks_file <- file.path(data_dir, "parks_properties.rds")
if (file.exists(parks_file)) {
parks_properties <- readRDS(parks_file)
} else {
parks_url <- "https://data.cityofnewyork.us/api/views/enfh-gkve/rows.csv?accessType=DOWNLOAD"
parks_properties <- read.csv(parks_url, stringsAsFactors = FALSE)
saveRDS(parks_properties, parks_file)
}
```
**Total parks loaded**: `r format(nrow(parks_properties), big.mark = ",")`
```{r process-parks}
#| label: process-parks
# Borough lookup table
borough_lookup <- c(
"M" = "Manhattan",
"B" = "Brooklyn",
"Q" = "Queens",
"X" = "Bronx",
"R" = "Staten Island"
)
# Process parks spatial data
parks_sf <- parks_properties %>%
filter(
!is.na(multipolygon),
multipolygon != "",
!is.na(BOROUGH),
BOROUGH != "",
BOROUGH %in% names(borough_lookup)
) %>%
st_as_sf(wkt = "multipolygon", crs = 4326) %>%
st_make_valid() %>%
mutate(
BOROUGH_FULL = recode(BOROUGH, !!!borough_lookup),
PROPNAME = SIGNNAME
) %>%
select(PROPNAME, BOROUGH = BOROUGH_FULL) %>%
st_transform(crs = 2263) # NYC State Plane (feet)
```
**Processed parks**: `r format(nrow(parks_sf), big.mark = ",")`
## U.S. Census Bureau Demographics Data
**Source**: American Community Survey 2022 5-Year Estimates
**API**: https://api.census.gov/data/2022/acs/acs5
Demographic data was obtained from the Census Bureau's ACS API, extracting six key variables (population, income, and four racial/ethnic categories) for New York State ZIP codes. Data was filtered to NYC using ZIP code ranges for the five boroughs, with records removed if population or income values were zero. The final dataset contains 177 NYC ZIP codes with complete demographic profiles.
### Variables Collected
- **B01003_001E**: Total Population
- **B19013_001E**: Median Household Income
- **B03002_003E**: White alone (Not Hispanic)
- **B03002_004E**: Black or African American alone (Not Hispanic)
- **B03002_006E**: Asian alone (Not Hispanic)
- **B03002_012E**: Hispanic or Latino
```{r load-census}
#| label: load-census
#| cache: true
demo_file <- file.path(data_dir, "nyc_census_demographics.rds")
if (file.exists(demo_file)) {
cat("Loading Census data from cache...\n")
nyc_demographics_raw <- readRDS(demo_file)
cat(sprintf("Loaded %d NYC ZIP codes\n", nrow(nyc_demographics_raw)))
} else {
cat("Downloading Census data from API...\n")
census_api_url <- paste0(
"https://api.census.gov/data/2022/acs/acs5?",
"get=NAME,B01003_001E,B19013_001E,B03002_003E,B03002_004E,B03002_006E,B03002_012E",
"&for=zip%20code%20tabulation%20area:*",
"&in=state:36"
)
tryCatch({
response <- GET(census_api_url, timeout(180))
if (status_code(response) != 200) {
stop(sprintf("Census API returned HTTP status code: %d", status_code(response)))
}
json_content <- content(response, as = "text", encoding = "UTF-8")
# Check for error messages in response
if (grepl("error|Error|ERROR", json_content)) {
cat("Census API returned an error. Response:\n")
cat(substr(json_content, 1, 500), "\n")
stop("Census API error - please run the R script first to cache the data")
}
census_raw <- fromJSON(json_content)
census_df <- as.data.frame(census_raw[-1, ], stringsAsFactors = FALSE)
colnames(census_df) <- census_raw[1, ]
census_df <- census_df %>%
mutate(
zip = `zip code tabulation area`,
total_pop = as.numeric(B01003_001E),
median_income = as.numeric(B19013_001E),
white = as.numeric(B03002_003E),
black = as.numeric(B03002_004E),
asian = as.numeric(B03002_006E),
hispanic = as.numeric(B03002_012E)
)
# Filter for NYC ZIP codes only
nyc_demographics_raw <- census_df %>%
mutate(zip_num = as.numeric(zip)) %>%
filter(
!is.na(zip_num),
(zip_num >= 10001 & zip_num <= 10282) |
(zip_num >= 10301 & zip_num <= 10314) |
(zip_num >= 10451 & zip_num <= 10475) |
(zip_num >= 11004 & zip_num <= 11109) |
(zip_num >= 11201 & zip_num <= 11256) |
(zip_num >= 11351 & zip_num <= 11697)
) %>%
filter(
!is.na(total_pop), total_pop > 0,
!is.na(median_income), median_income > 0
) %>%
mutate(
pct_white = (white / total_pop) * 100,
pct_black = (black / total_pop) * 100,
pct_asian = (asian / total_pop) * 100,
pct_hispanic = (hispanic / total_pop) * 100,
borough = case_when(
zip_num >= 10001 & zip_num <= 10282 ~ "Manhattan",
zip_num >= 10301 & zip_num <= 10314 ~ "Staten Island",
zip_num >= 10451 & zip_num <= 10475 ~ "Bronx",
zip_num >= 11201 & zip_num <= 11256 ~ "Brooklyn",
zip_num >= 11004 & zip_num <= 11109 ~ "Queens",
zip_num >= 11351 & zip_num <= 11697 ~ "Queens",
TRUE ~ NA_character_
)
) %>%
filter(!is.na(borough)) %>%
select(zip, borough, total_pop, median_income,
pct_white, pct_black, pct_asian, pct_hispanic)
saveRDS(nyc_demographics_raw, demo_file)
cat(sprintf("Downloaded and cached %d NYC ZIP codes\n", nrow(nyc_demographics_raw)))
}, error = function(e) {
cat("\n=======================================================================\n")
cat("ERROR: Census API request failed\n")
cat("=======================================================================\n")
cat("Error message:", e$message, "\n\n")
cat("SOLUTION: Please run the R script first to cache the data:\n")
cat(" source('final_individual_report.R')\n\n")
cat("This will download and cache all data, then the QMD will work.\n")
cat("=======================================================================\n\n")
stop("Census data not available. Run R script first to cache data.")
})
}
```
**Total NYC ZIP codes**: `r nrow(nyc_demographics_raw)`
## NYC ZIP Code Boundaries
**Source**: NYC Open Data - Modified ZIP Code Tabulation Areas (MODZCTA)
**URL**: https://data.cityofnewyork.us/resource/pri4-ifjk.geojson
The MODZCTA dataset provides GeoJSON polygon geometries for NYC residential ZIP code areas, modified to better reflect actual neighborhoods by excluding non-residential areas. Data was imported using `st_read()` and transformed from WGS84 to NYC State Plane (EPSG:2263) for accurate distance measurements.
```{r load-boundaries}
#| label: load-boundaries
zip_boundaries_file <- file.path(data_dir, "zip_boundaries.rds")
if (file.exists(zip_boundaries_file)) {
zip_boundaries <- readRDS(zip_boundaries_file)
} else {
zip_url <- "https://data.cityofnewyork.us/resource/pri4-ifjk.geojson"
zip_boundaries <- st_read(zip_url, quiet = TRUE)
zip_boundaries <- st_transform(zip_boundaries, crs = 2263)
saveRDS(zip_boundaries, zip_boundaries_file)
}
```
**Total ZIP boundaries**: `r nrow(zip_boundaries)`
## Merging and Final Data Preparation
Census demographics were joined to MODZCTA spatial boundaries using ZIP codes as the join key. Both datasets required conversion to character format for successful matching. Income categories were created using thresholds ($40k, $75k, $125k), and majority demographic groups were identified when any racial/ethnic category exceeded 50% of population. The final dataset contains 177 ZIP code areas representing approximately 8.3 million NYC residents.
```{r merge-data}
#| label: merge-data
zip_col_name <- intersect(
names(zip_boundaries),
c("modzcta", "MODZCTA", "ZIPCODE", "zipcode", "zip", "ZIP", "zcta", "ZCTA")
)[1]
zip_boundaries_clean <- zip_boundaries %>%
mutate(zip_from_boundary = as.character(.data[[zip_col_name]]))
nyc_demographics_sf <- zip_boundaries_clean %>%
left_join(
nyc_demographics_raw %>%
mutate(zip = as.character(zip)) %>%
select(zip, borough, total_pop, median_income,
pct_white, pct_black, pct_hispanic, pct_asian),
by = c("zip_from_boundary" = "zip"),
relationship = "many-to-one"
) %>%
filter(
!is.na(total_pop), total_pop > 0,
!is.na(median_income), median_income > 0
)
# Create categorical variables
nyc_demographics_sf <- nyc_demographics_sf %>%
mutate(
zip = zip_from_boundary,
income_category = case_when(
median_income < 40000 ~ "Low Income (<$40k)",
median_income < 75000 ~ "Lower-Middle ($40k-$75k)",
median_income < 125000 ~ "Upper-Middle ($75k-$125k)",
median_income >= 125000 ~ "High Income (>$125k)",
TRUE ~ "Unknown"
),
income_category = factor(income_category, levels = c(
"Low Income (<$40k)",
"Lower-Middle ($40k-$75k)",
"Upper-Middle ($75k-$125k)",
"High Income (>$125k)"
)),
majority_group = case_when(
pct_hispanic > 50 ~ "Majority Hispanic",
pct_white > 50 ~ "Majority White",
pct_black > 50 ~ "Majority Black",
pct_asian > 50 ~ "Majority Asian",
TRUE ~ "No Majority"
)
)
```
**Successfully merged**: `r nrow(nyc_demographics_sf)` ZIP areas with complete data
---
# Distance Calculation Methodology {#sec-methodology}
## Approach: Euclidean (Straight-Line) Distance
We calculate the shortest straight-line distance from each ZIP code centroid to the nearest park boundary. This method provides a consistent, reproducible measure of park proximity.
### Step 1: Calculate ZIP Code Centroids
```{r calculate-centroids}
#| label: calculate-centroids
zip_centroids <- st_centroid(nyc_demographics_sf)
```
- Computed centroids for **`r nrow(zip_centroids)` ZIP code areas**
- Centroids represent the geographic center of each ZIP area
### Step 2: Find Nearest Park
```{r find-nearest}
#| label: find-nearest
nearest_park_idx <- st_nearest_feature(zip_centroids, parks_sf)
```
- Used spatial indexing for efficient nearest-neighbor search
- Each ZIP code matched to its geographically closest park
### Step 3: Calculate Distances
```{r calculate-distances}
#| label: calculate-distances
distances <- st_distance(zip_centroids, parks_sf[nearest_park_idx, ], by_element = TRUE)
distances_mi <- as.numeric(distances) / 5280 # Convert feet to miles
```
- Distances calculated in feet using NYC State Plane projection
- Converted to miles for interpretability (1 mile = 5,280 feet)
- Distance range: **`r sprintf("%.3f - %.3f miles", min(distances_mi), max(distances_mi))`**
### Methodological Limitation
::: {.callout-note icon=false}
## Important Note
Euclidean distance is a proxy for actual walking distance. True walking paths would need to account for street networks, pedestrian infrastructure, and barriers (highways, rivers). Research suggests straight-line distance underestimates actual walking distance by approximately 20-30%.
:::
### Benchmark: 10-Minute Walk Standard
The Trust for Public Land defines park accessibility as being within a **10-minute walk (approximately 0.5 miles)**. We use this as our equity benchmark.
```{r prepare-results}
#| label: prepare-results
distance_results <- nyc_demographics_sf %>%
st_drop_geometry() %>%
mutate(
distance_to_park_ft = as.numeric(distances),
distance_to_park_mi = distances_mi,
nearest_park_name = parks_sf$PROPNAME[nearest_park_idx],
nearest_park_borough = parks_sf$BOROUGH[nearest_park_idx],
within_10min_walk = distance_to_park_mi <= 0.5
)
pct_accessible <- mean(distance_results$within_10min_walk, na.rm = TRUE) * 100
```
**Preliminary Finding**: `r sprintf("%.1f%%", pct_accessible)` of NYC ZIP codes meet the 10-minute walk standard.
---
# Data Analysis and Visualizations {#sec-analysis}
```{r theme-setup}
#| label: theme-setup
# Professional theme for all plots
theme_professional <- function() {
theme_minimal(base_size = 13, base_family = "sans") +
theme(
plot.title = element_text(face = "bold", size = 16, hjust = 0, margin = margin(b = 8)),
plot.subtitle = element_text(size = 11, color = "gray40", hjust = 0, margin = margin(b = 15)),
plot.caption = element_text(size = 9, color = "gray50", hjust = 1, margin = margin(t = 15)),
panel.grid.major = element_line(color = "gray90", linewidth = 0.3),
panel.grid.minor = element_blank(),
panel.background = element_rect(fill = "white", color = NA),
plot.background = element_rect(fill = "white", color = NA),
axis.title = element_text(size = 11, face = "bold", color = "gray30"),
axis.text = element_text(size = 10, color = "gray40"),
axis.line = element_line(color = "gray70", linewidth = 0.3),
axis.ticks = element_line(color = "gray70", linewidth = 0.3),
legend.position = "right",
legend.title = element_text(size = 10, face = "bold"),
legend.text = element_text(size = 9),
legend.background = element_rect(fill = "white", color = "gray80", linewidth = 0.3),
plot.margin = margin(20, 20, 20, 20)
)
}
theme_set(theme_professional())
# Custom color palettes
income_colors <- c(
"Low Income (<$40k)" = "#d73027",
"Lower-Middle ($40k-$75k)" = "#fc8d59",
"Upper-Middle ($75k-$125k)" = "#91bfdb",
"High Income (>$125k)" = "#4575b4"
)
demographic_colors <- c(
"Majority Hispanic" = "#e78ac3",
"Majority White" = "#8da0cb",
"Majority Black" = "#fc8d62",
"Majority Asian" = "#66c2a5",
"No Majority" = "#a6d854"
)
```
## Visualization 1: Interactive Spatial Distribution Map
```{r interactive-map}
#| label: fig-interactive-map
#| fig-cap: "Interactive map showing walking distance to nearest park by NYC ZIP code. Click on areas for detailed demographic information."
#| fig-width: 10
#| fig-height: 8
# Prepare map data
map_data <- nyc_demographics_sf %>%
mutate(
distance_to_park_mi = distance_results$distance_to_park_mi,
distance_to_park_ft = distance_results$distance_to_park_ft,
nearest_park_name = distance_results$nearest_park_name,
nearest_park_borough = distance_results$nearest_park_borough,
within_10min_walk = distance_results$within_10min_walk
)
# Transform to WGS84 for Leaflet
map_data_wgs84 <- st_transform(map_data, crs = 4326)
# Create color palette
pal_distance <- colorNumeric(
palette = c("#065f46", "#10b981", "#fbbf24", "#f59e0b", "#dc2626"),
domain = map_data_wgs84$distance_to_park_mi,
na.color = "#808080"
)
# Create popup content
popup_content <- paste0(
"<div style='font-family: Arial; font-size: 12px; padding: 5px;'>",
"<b style='font-size: 14px;'>ZIP Code: ", map_data_wgs84$zip, "</b><br><br>",
"<b>Park Access Metrics:</b><br>",
"Distance to Park: <b>", round(map_data_wgs84$distance_to_park_mi, 3), " miles</b><br>",
"Nearest Park: ", map_data_wgs84$nearest_park_name, "<br>",
"10-Min Walk: <b>", ifelse(map_data_wgs84$within_10min_walk,
"<span style='color: green;'>✓ Yes</span>",
"<span style='color: red;'>✗ No</span>"), "</b><br><br>",
"<b>Demographics:</b><br>",
"Borough: ", map_data_wgs84$borough, "<br>",
"Population: ", format(round(map_data_wgs84$total_pop), big.mark = ","), "<br>",
"Median Income: $", format(round(map_data_wgs84$median_income), big.mark = ","), "<br>",
"Income Category: ", map_data_wgs84$income_category, "<br><br>",
"<b>Racial/Ethnic Composition:</b><br>",
"White: ", round(map_data_wgs84$pct_white, 1), "%<br>",
"Black: ", round(map_data_wgs84$pct_black, 1), "%<br>",
"Hispanic: ", round(map_data_wgs84$pct_hispanic, 1), "%<br>",
"Asian: ", round(map_data_wgs84$pct_asian, 1), "%<br>",
"</div>"
)
# Create interactive map
leaflet(map_data_wgs84, width = "100%", height = "600px") %>%
addProviderTiles(providers$CartoDB.Positron) %>%
addPolygons(
fillColor = ~pal_distance(distance_to_park_mi),
fillOpacity = 0.7,
color = "white",
weight = 1,
popup = popup_content,
highlight = highlightOptions(
weight = 3,
color = "#1e3a8a",
fillOpacity = 0.9,
bringToFront = TRUE
),
label = ~paste0("ZIP ", zip, ": ", round(distance_to_park_mi, 3), " miles"),
labelOptions = labelOptions(
style = list("font-weight" = "normal", padding = "3px 8px"),
textsize = "12px",
direction = "auto"
)
) %>%
addLegend(
position = "bottomright",
pal = pal_distance,
values = ~distance_to_park_mi,
title = "Distance to<br>Nearest Park<br>(miles)",
opacity = 0.8,
labFormat = labelFormat(digits = 2)
) %>%
setView(lng = -73.935, lat = 40.730, zoom = 10)
```
### Visualization Design
This **interactive map** employs a **5-color gradient scheme** to represent walking distance to the nearest park across NYC ZIP codes. The color palette ranges from deep green (excellent access, <0.25 miles) through yellow/orange (moderate access) to dark red (poor access, >1.0 mile).
### Color Scheme Rationale
- **Dark Green (#065f46)**: < 0.25 mi - Excellent access
- **Green (#10b981)**: 0.25-0.5 mi - Good access (within 10-min walk)
- **Yellow (#fbbf24)**: 0.5-0.7 mi - Moderate access
- **Orange (#f59e0b)**: 0.7-1.0 mi - Limited access
- **Red (#dc2626)**: > 1.0 mi - Poor access
### Distance Calculation Method
As described in @sec-methodology, distances are calculated as:
1. ZIP code centroid coordinates (geographic center of area)
2. Nearest park boundary (using spatial indexing)
3. Euclidean distance in NYC State Plane projection (feet)
4. Conversion to miles for interpretability
This method provides consistent, reproducible measurements but underestimates true walking distance by ~20-30% due to street network constraints.
### Key Spatial Patterns
```{r borough-stats}
#| label: tbl-borough-stats
#| tbl-cap: "Park access statistics by borough"
borough_stats <- distance_results %>%
group_by(borough) %>%
summarise(
`ZIP Codes` = n(),
`Avg Distance (mi)` = round(mean(distance_to_park_mi, na.rm = TRUE), 3),
`Median Distance (mi)` = round(median(distance_to_park_mi, na.rm = TRUE), 3),
`% Within 10-Min Walk` = round(mean(within_10min_walk, na.rm = TRUE) * 100, 1),
.groups = "drop"
) %>%
arrange(desc(`Avg Distance (mi)`))
knitr::kable(borough_stats)
```
---
## Visualization 2: Income-Distance Relationship Analysis
```{r correlation-data-prep}
#| label: correlation-data-prep
cor_data <- distance_results %>%
filter(complete.cases(distance_to_park_mi, median_income))
cor_test <- cor.test(cor_data$distance_to_park_mi, cor_data$median_income)
# Linear regression model
lm_model <- lm(distance_to_park_mi ~ median_income, data = cor_data)
lm_summary <- summary(lm_model)
```
```{r scatter-plot}
#| label: fig-scatter
#| fig-cap: "Relationship between neighborhood median household income and walking distance to nearest park. Each point represents a ZIP code, sized by population. Red line shows linear regression with 95% confidence interval."
#| fig-width: 10
#| fig-height: 7
ggplot(cor_data, aes(x = median_income, y = distance_to_park_mi)) +
geom_point(
aes(color = income_category, size = total_pop),
alpha = 0.6
) +
geom_smooth(
method = "lm",
se = TRUE,
color = "#dc2626",
fill = "#fecaca",
linetype = "solid",
linewidth = 1.2
) +
scale_x_continuous(
labels = dollar_format(),
breaks = seq(0, max(cor_data$median_income, na.rm = TRUE), 25000)
) +
scale_y_continuous(
breaks = seq(0, max(cor_data$distance_to_park_mi, na.rm = TRUE), 0.2)
) +
scale_color_manual(
values = income_colors,
name = "Income Category"
) +
scale_size_continuous(
range = c(2, 10),
labels = comma_format()
) +
labs(
title = "Relationship Between Neighborhood Income and Park Access",
subtitle = sprintf(
"Pearson r = %.3f (p %s) | R² = %.3f | n = %d ZIP codes",
cor_test$estimate,
ifelse(cor_test$p.value < 0.001, "< 0.001", sprintf("= %.3f", cor_test$p.value)),
lm_summary$r.squared,
nrow(cor_data)
),
x = "Median Household Income",
y = "Walking Distance to Nearest Park (miles)",
caption = paste0(
"Linear regression line with 95% confidence interval (shaded area)\n",
"Point size represents ZIP code population | Data: NYC Parks Properties & Census ACS 2022"
),
size = "Population"
) +
theme_professional() +
theme(
legend.position = "right",
legend.box = "vertical"
)
```
### Statistical Analysis
We examine the relationship between neighborhood median household income and walking distance to the nearest park using **Pearson correlation** and **linear regression modeling**.
### Correlation Analysis Results
- **Pearson's r**: `r sprintf("%.4f", cor_test$estimate)`
- **p-value**: `r sprintf("%.6f", cor_test$p.value)`
- **95% CI**: [`r sprintf("%.4f", cor_test$conf.int[1])`, `r sprintf("%.4f", cor_test$conf.int[2])`]
- **Sample size**: `r nrow(cor_data)` ZIP codes
### Interpretation
```{r interpretation}
#| label: interpretation
#| echo: false
if (cor_test$p.value < 0.001) {
sig_level <- "highly significant (p < 0.001)"
} else if (cor_test$p.value < 0.01) {
sig_level <- "very significant (p < 0.01)"
} else if (cor_test$p.value < 0.05) {
sig_level <- "statistically significant (p < 0.05)"
} else {
sig_level <- "not statistically significant (p ≥ 0.05)"
}
if (abs(cor_test$estimate) < 0.3) {
effect_size <- "weak"
} else if (abs(cor_test$estimate) < 0.5) {
effect_size <- "moderate"
} else {
effect_size <- "strong"
}
```
The correlation is **`r sig_level`**.
`r if (cor_test$estimate < 0) paste0("**Negative correlation** indicates that as neighborhood income increases, walking distance to parks DECREASES. This suggests higher-income neighborhoods have better park access - evidence of environmental inequity.")`
`r if (cor_test$estimate >= 0) paste0("**Positive correlation** indicates that as neighborhood income increases, walking distance to parks INCREASES. This would suggest lower-income neighborhoods have better access - contrary to typical environmental justice patterns.")`
**Effect size**: `r effect_size` (|r| = `r sprintf("%.3f", abs(cor_test$estimate))`)
### Linear Regression Model
- **Model**: Distance = `r sprintf("%.6f", coef(lm_model)[1])` + `r sprintf("%.9f", coef(lm_model)[2])` × Income
- **R-squared**: `r sprintf("%.4f", lm_summary$r.squared)` (`r sprintf("%.1f%%", lm_summary$r.squared * 100)` of variance explained)
- **F-statistic**: `r sprintf("%.2f", lm_summary$fstatistic[1])` (p = `r sprintf("%.6f", pf(lm_summary$fstatistic[1], lm_summary$fstatistic[2], lm_summary$fstatistic[3], lower.tail = FALSE))`)
### Practical Significance
```{r practical-sig}
#| label: practical-sig
#| echo: false
income_10k_increase <- coef(lm_model)[2] * 10000
```
For every **$10,000 increase** in median income, walking distance changes by **`r sprintf("%.4f", income_10k_increase)` miles** (`r sprintf("%.0f", income_10k_increase * 5280)` feet).
`r if (income_10k_increase < 0) paste0("This **negative relationship confirms environmental inequity**: wealthier neighborhoods systematically enjoy closer proximity to parks.")`
### Statistical Validity Assessment
**Regression diagnostics:**
- Linear relationship: Visually confirmed in scatter plot
- Independence: ZIP codes are independent geographic units
- Sample size: Adequate for reliable inference (n > 30)
::: {.callout-note icon=false}
## Conclusion
`r if (cor_test$p.value < 0.05 && abs(cor_test$estimate) > 0.1) paste0("There is **statistically significant evidence** of a relationship between neighborhood income and park access. This finding supports the hypothesis that park accessibility is not equitably distributed across NYC neighborhoods.")`
`r if (cor_test$p.value >= 0.05 || abs(cor_test$estimate) <= 0.1) paste0("The relationship between income and park access is not statistically significant or is very weak, suggesting income may not be a primary determinant of park proximity in NYC.")`
:::
---
## Visualization 3: Park Access by Income Category
```{r income-data-prep}
#| label: income-data-prep
income_analysis <- distance_results %>%
filter(!is.na(income_category)) %>%
group_by(income_category) %>%
summarise(
n_areas = n(),
total_pop = sum(total_pop, na.rm = TRUE),
avg_distance_mi = mean(distance_to_park_mi, na.rm = TRUE),
median_distance_mi = median(distance_to_park_mi, na.rm = TRUE),
sd_distance_mi = sd(distance_to_park_mi, na.rm = TRUE),
se_distance_mi = sd_distance_mi / sqrt(n_areas),
pct_accessible = mean(within_10min_walk, na.rm = TRUE) * 100,
.groups = "drop"
)
```
```{r income-plot}
#| label: fig-income
#| fig-cap: "Average walking distance to parks by neighborhood income level. Error bars represent standard error. Dashed line shows 10-minute walk threshold (0.5 miles)."
#| fig-width: 10
#| fig-height: 7
ggplot(income_analysis,
aes(x = income_category, y = avg_distance_mi, fill = income_category)) +
geom_col(width = 0.7, alpha = 0.9) +
geom_errorbar(
aes(ymin = avg_distance_mi - se_distance_mi,
ymax = avg_distance_mi + se_distance_mi),
width = 0.25,
linewidth = 0.6,
color = "gray30"
) +
geom_hline(
yintercept = 0.5,
linetype = "dashed",
color = "#059669",
linewidth = 1
) +
geom_text(
aes(label = sprintf("%.3f mi\n(%d ZIPs)", avg_distance_mi, n_areas)),
vjust = -0.5,
size = 3.5,
fontface = "bold",
color = "gray30"
) +
annotate(
"text",
x = 4.5,
y = 0.5,
label = "10-min walk threshold",
vjust = -0.5,
hjust = 1,
size = 3,
color = "#059669",
fontface = "italic"
) +
scale_fill_manual(
values = income_colors,
guide = "none"
) +
scale_y_continuous(
expand = expansion(mult = c(0, 0.15)),
breaks = seq(0, 1, 0.1)
) +
labs(
title = "Average Walking Distance to Parks by Neighborhood Income Level",
subtitle = "Error bars represent standard error of the mean",
x = NULL,
y = "Average Distance to Nearest Park (miles)",
caption = "Dashed line indicates 10-minute walk threshold (0.5 miles)"
) +
theme_professional() +
theme(
axis.text.x = element_text(angle = 0, hjust = 0.5, size = 10, face = "bold"),
panel.grid.major.x = element_blank()
)
```
```{r income-table}
#| label: tbl-income
#| tbl-cap: "Park access metrics by income category"
knitr::kable(income_analysis, digits = 3)
```
### Statistical Comparison of Income Groups
```{r anova}
#| label: anova-test
# ANOVA test
anova_test <- aov(distance_to_park_mi ~ income_category, data = distance_results)
anova_summary <- summary(anova_test)
```
**One-Way ANOVA Results:**
- F-statistic: `r sprintf("%.3f", anova_summary[[1]]$'F value'[1])`
- p-value: `r sprintf("%.6f", anova_summary[[1]]$'Pr(>F)'[1])`
`r if (anova_summary[[1]]$'Pr(>F)'[1] < 0.05) paste0("**Significant differences exist** between income groups (p < 0.05).")`
---
## Visualization 4: Park Access by Racial/Ethnic Composition
```{r race-analysis}
#| label: race-analysis
race_analysis <- distance_results %>%
group_by(majority_group) %>%
summarise(
n_areas = n(),
total_pop = sum(total_pop, na.rm = TRUE),
avg_distance_mi = mean(distance_to_park_mi, na.rm = TRUE),
median_distance_mi = median(distance_to_park_mi, na.rm = TRUE),
sd_distance_mi = sd(distance_to_park_mi, na.rm = TRUE),
se_distance_mi = sd_distance_mi / sqrt(n_areas),
pct_accessible = mean(within_10min_walk, na.rm = TRUE) * 100,
avg_income = mean(median_income, na.rm = TRUE),
.groups = "drop"
) %>%
arrange(desc(avg_distance_mi))
```
```{r race-plot}
#| label: fig-race
#| fig-cap: "Average walking distance to parks by neighborhood racial/ethnic composition. Neighborhoods grouped by majority (>50%) demographic group."
#| fig-width: 10
#| fig-height: 7
ggplot(race_analysis,
aes(x = reorder(majority_group, -avg_distance_mi),
y = avg_distance_mi,
fill = majority_group)) +
geom_col(width = 0.7, alpha = 0.9) +
geom_errorbar(
aes(ymin = avg_distance_mi - se_distance_mi,
ymax = avg_distance_mi + se_distance_mi),
width = 0.25,
linewidth = 0.6,
color = "gray30"
) +
geom_hline(
yintercept = 0.5,
linetype = "dashed",
color = "#059669",
linewidth = 1
) +
geom_text(
aes(label = sprintf("%.3f mi", avg_distance_mi)),
vjust = -2.5,
size = 3.5,
fontface = "bold",
color = "gray30"
) +
geom_text(
aes(label = sprintf("n=%d", n_areas)),
vjust = 3,
size = 3,
color = "white",
fontface = "bold"
) +
scale_fill_manual(
values = demographic_colors,
guide = "none"
) +
scale_y_continuous(
expand = expansion(mult = c(0, 0.3)),
breaks = seq(0, 1, 0.1)
) +
labs(
title = "Average Walking Distance to Parks by Neighborhood Demographics",
subtitle = "Grouped by majority racial/ethnic composition (>50% of population)",
x = NULL,
y = "Average Distance to Nearest Park (miles)",
caption = "Error bars represent standard error | Dashed line indicates 10-minute walk threshold"
) +
theme_professional() +
theme(
axis.text.x = element_text(angle = 20, hjust = 1, size = 10, face = "bold"),
panel.grid.major.x = element_blank()
)
```
```{r race-table}
#| label: tbl-race
#| tbl-cap: "Park access metrics by racial/ethnic composition"
knitr::kable(race_analysis, digits = 3)
```
---
# Key Findings and Insights {#sec-findings}
## Answering the Specific Question
**How does the average walking distance to the nearest park vary across the city and neighborhood demographic composition?**
### Finding 1: Overall Park Accessibility
```{r overall-stats}
#| label: overall-stats
#| echo: false
mean_dist <- mean(distance_results$distance_to_park_mi)
median_dist <- median(distance_results$distance_to_park_mi)
total_pop <- sum(distance_results$total_pop)
```
- **Mean distance**: `r sprintf("%.3f", mean_dist)` miles
- **Median distance**: `r sprintf("%.3f", median_dist)` miles
- **ZIP codes within 10-min walk**: `r sprintf("%.1f%%", pct_accessible)`
- **ZIP codes analyzed**: `r nrow(distance_results)`
- **Total population represented**: `r format(total_pop, big.mark = ",")`
::: {.callout-note icon=false}
## Interpretation
`r if (pct_accessible < 50) paste0("**Less than half** of NYC ZIP codes meet the 10-minute walk standard, indicating widespread accessibility challenges.")`
`r if (pct_accessible >= 50 && pct_accessible < 75) paste0("While a **majority** of ZIP codes meet the 10-minute walk standard, significant gaps remain for many neighborhoods.")`
`r if (pct_accessible >= 75) paste0("**Most** NYC ZIP codes meet the 10-minute walk standard, though disparities exist for specific communities.")`
:::
### Finding 2: Income-Based Disparities
```{r income-disparity}
#| label: income-disparity
#| echo: false
low_income <- income_analysis %>% filter(income_category == "Low Income (<$40k)")
high_income <- income_analysis %>% filter(income_category == "High Income (>$125k)")
if (nrow(low_income) > 0 && nrow(high_income) > 0) {
disparity_ratio <- low_income$avg_distance_mi / high_income$avg_distance_mi
disparity_pct <- ((low_income$avg_distance_mi - high_income$avg_distance_mi) /
high_income$avg_distance_mi) * 100
}
```
`r if (nrow(low_income) > 0 && nrow(high_income) > 0) sprintf("- **Low income areas**: %.3f miles average\n- **High income areas**: %.3f miles average\n- **Disparity ratio**: %.2fx\n- **Percent difference**: %.1f%%", low_income$avg_distance_mi, high_income$avg_distance_mi, disparity_ratio, abs(disparity_pct))`
::: {.callout-note icon=false}
## Interpretation
`r if (exists("disparity_ratio") && disparity_ratio > 1.2) paste0("**Low-income neighborhoods face significantly worse park access** than high-income areas, demonstrating clear environmental inequity.")`
`r if (exists("disparity_ratio") && disparity_ratio < 0.8) paste0("Surprisingly, **low-income neighborhoods have better park access** than high-income areas, contradicting typical environmental justice patterns.")`
`r if (exists("disparity_ratio") && disparity_ratio >= 0.8 && disparity_ratio <= 1.2) paste0("Park access is **relatively equitable** across income levels, with minimal differences between low and high-income areas.")`
:::
### Finding 3: Racial/Ethnic Disparities
```{r race-disparity}
#| label: race-disparity
#| echo: false
if (nrow(race_analysis) > 0) {
worst_access <- race_analysis %>% slice_max(avg_distance_mi, n = 1)
best_access <- race_analysis %>% slice_min(avg_distance_mi, n = 1)
}
```
`r if (nrow(race_analysis) > 0) sprintf("- **Worst access**: %s (%.3f miles)\n- **Best access**: %s (%.3f miles)\n- **Disparity ratio**: %.2fx", worst_access$majority_group, worst_access$avg_distance_mi, best_access$majority_group, best_access$avg_distance_mi, worst_access$avg_distance_mi / best_access$avg_distance_mi)`
Racial and ethnic composition shows variation in park access, suggesting that environmental inequity has demographic dimensions beyond income alone.
### Finding 4: Statistical Significance
- **Income-distance correlation**: r = `r sprintf("%.3f", cor_test$estimate)` (p = `r sprintf("%.4f", cor_test$p.value)`)
- **Effect size**: `r effect_size`
- **R-squared**: `r sprintf("%.3f", lm_summary$r.squared)` (`r sprintf("%.1f%%", lm_summary$r.squared * 100)` variance explained)
`r if (cor_test$p.value < 0.05) paste0("Statistical tests **confirm a significant relationship** between neighborhood characteristics and park access, providing evidence-based support for park equity interventions.")`
## Contribution to Overarching Question
**How equitable is access to public park space across different neighborhoods in New York City?**
Our proximity analysis reveals that **park equity cannot be assessed by quantity alone**. Key insights include:
### 1. Spatial Inequality
Geographic distance creates access barriers even where parks exist. Our findings show systematic variation in walking distances across demographic groups.
### 2. Intersectional Disparities
Both income and race/ethnicity correlate with park access, suggesting **multiple dimensions of inequity** that require targeted interventions.
### 3. Policy Implications
Meeting the 10-minute walk standard requires both:
- New park development in underserved areas
- Improved pedestrian infrastructure to reduce effective walking distances
### 4. Integration with Team Findings
When combined with **quantity** (acreage per capita) and **quality** (amenity assessments), proximity metrics complete a comprehensive equity framework. Some neighborhoods may score well on one dimension but poorly on others, requiring multifaceted interventions.
---
# Conclusions {#sec-conclusions}
This spatial analysis provides quantitative evidence that park access in New York City exhibits significant spatial and demographic disparities. Our investigation reveals that only `r sprintf("%.1f%%", pct_accessible)` of NYC ZIP codes meet the 10-minute walk standard for park accessibility, with clear patterns emerging across income and racial/ethnic lines. `r if (cor_test$estimate < 0 && cor_test$p.value < 0.05) paste0("The statistically significant negative correlation between neighborhood income and park distance (r = ", sprintf("%.3f", cor_test$estimate), ", p < 0.05) demonstrates that higher-income neighborhoods systematically enjoy closer proximity to parks, while lower-income communities face longer walking distances - a clear manifestation of environmental inequity. Similarly, racial and ethnic composition correlates with access disparities, with ", if(nrow(race_analysis) > 0) sprintf("%s neighborhoods experiencing the poorest access at %.3f miles average distance", race_analysis %>% slice_max(avg_distance_mi, n=1) %>% pull(majority_group), race_analysis %>% slice_max(avg_distance_mi, n=1) %>% pull(avg_distance_mi)), " compared to better-served areas.")`
Combined with team members' research on park quantity and quality, our proximity analysis completes a comprehensive view of environmental justice issues in urban green space, revealing that some neighborhoods face compounding disadvantages - greater distances to parks, less total park space per capita, and lower-quality amenities.
The path toward equity requires multi-dimensional interventions: building new parks in underserved areas identified through our spatial analysis, improving pedestrian infrastructure to reduce effective walking distances, ensuring park quality meets community needs, and maintaining accountability through ongoing measurement and reporting. By addressing these spatial, quantitative, and qualitative dimensions together, New York City can move toward ensuring all residents - regardless of income, race, or neighborhood - enjoy equitable access to high-quality public park space, recognizing that true environmental justice demands both proximity and quality in park provision.
---
# References {#sec-references}
## Data Sources
[1] **NYC Parks Properties**
Source: NYC Open Data
URL: https://data.cityofnewyork.us/api/views/enfh-gkve/rows.csv
Accessed: December 2025
Records: `r format(nrow(parks_properties), big.mark = ",")` parks
[2] **U.S. Census Bureau - American Community Survey 2022 (5-Year Estimates)**
Source: U.S. Census Bureau API
URL: https://api.census.gov/data/2022/acs/acs5
Accessed: December 2025
Variables: B01003_001E, B19013_001E, B03002_003E, B03002_004E, B03002_006E, B03002_012E
[3] **NYC Modified ZIP Code Tabulation Areas (MODZCTA)**
Source: NYC Open Data
URL: https://data.cityofnewyork.us/resource/pri4-ifjk.geojson
Accessed: December 2025
Records: `r format(nrow(zip_boundaries), big.mark = ",")` ZIP boundaries
## Methodology References
[4] Trust for Public Land (2023). *ParkScore Index Methodology*. 10-minute walk standard (0.5 miles) benchmark.
[5] Sister, C., Wolch, J., & Wilson, J. (2010). Got green? Addressing environmental justice in park provision. *GeoJournal*, 75(3), 229-248.
[6] Rigolon, A. (2016). A complex landscape of inequity in access to urban parks: A literature review. *Landscape and Urban Planning*, 153, 160-169.
---
